Chapter 16

The aim of these exercises is to continue making sure you’re comfortable with handling multivariate data. In this chapter’s, you’ll focus on analyses based on dissimilarities/distances, including fitting linear models to these kinds of response variables.

For each of the examples below, you should follow the sequence we’ve used previously, as far as it’s sensible:

1. What is the biological question?

2. Is the predictor continuous or categorical?

3. Write out the linear model corresponding to this question.

4. What distribution do you expect the response variable to follow?

5. What are the assumptions behind the statistical model you’ll fit?

   1. Are those assumptions satisfied?

6. Fit the model

   1. How will you assess whether the model fits well?

   2. Can you detect an effect of the predictors?

   3. How do you measure the effect?

7. What do you conclude (including any cautions)

---

A

Dixon et al. (2018) focused on assemblages of reptiles in forests and woodlands of southeastern Australia, with a particular interest in relationships between these assemblages and the fire history of the landscape. They identified 81 sites that varied in time since fire, from 6 months to >96 y. Rather than a continuum, three categories of time since fire were used, 0.5-2y, 6-12y, and >96y. Sites were also classified according to habitat.

Reptile assemblages were sampled with a range of methods, including visual surveys and camera traps, to give counts of 20 reptiles, whose abundance ranged from 0-1 to 0-126, depending on species.

Data are available from dryad, as Rep_abund.csv. You’ll want to focus on the columns with reptile numbers, and the one with the fire history (tsf). For convenience, we’ve extracted those data as dixonbiota.csv

Start by having a quick look at the data.

df <- read_csv("../data/dixonbiota.csv")

Rows: 79 Columns: 22── Column specification ────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): site, tsf
dbl (20): adup, aplat, amur, amac, aram, dcor, ecun, esax, eul, htal, ldel, lgui, lwhi, ppor...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(df, 10)

The file is pretty straightforward, with tsf being the fire category, and columns 3-22 each representing one reptile taxon.

Outline how you’d assess whether the reptile assemblage (the combination of species present and their abundances) is related to fire history, using a dissimilarity-based approach.

You should think about how you’ll deal with the data:

• Will you need a transformation, standardization, etc.?

• What are the consequences of any decisions you make for the interpretation of your analysis?

• What measure(s) of dissimilarity are appropriate?

Run the analysis and provide your interpretation

MDS on reptile abundances

Start with dissimilarities

Students should have some discussion of raw vs transformed (e.g. 4th root), standardization.

Fourth root is common. Using it downweights the influence of abundant species, while raw allows those species to have a strong influence. There are questions about which one of these reflects ecological function better

Data are all abundances

df1 <- df[,-(1:2)]   #Stripped back file without labels, leaving only 
df1s <- wisconsin(df1)   #Common to do some standardisation
df1s.bc <- vegdist(df1s,'bray')   #With standardisation (wisconsin)

Visualize results

Need to decide whether to use 2 or 3 dimensions

#Run 3d first
df1s.mds3 <- metaMDS(df1s.bc,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.184364 
Run 1 stress 0.1835256 
... New best solution
... Procrustes: rmse 0.02423641  max resid 0.1643187 
Run 2 stress 0.1835251 
... New best solution
... Procrustes: rmse 0.000295326  max resid 0.001495451 
... Similar to previous best
Run 3 stress 0.1880849 
Run 4 stress 0.1939628 
Run 5 stress 0.1873871 
Run 6 stress 0.1835251 
... Procrustes: rmse 0.0004838384  max resid 0.002516706 
... Similar to previous best
Run 7 stress 0.1861646 
Run 8 stress 0.1835248 
... New best solution
... Procrustes: rmse 0.0003669419  max resid 0.002417968 
... Similar to previous best
Run 9 stress 0.1836813 
... Procrustes: rmse 0.007408813  max resid 0.04927581 
Run 10 stress 0.1882458 
Run 11 stress 0.1835256 
... Procrustes: rmse 0.0005098868  max resid 0.003722284 
... Similar to previous best
Run 12 stress 0.1843908 
Run 13 stress 0.1874228 
Run 14 stress 0.1835261 
... Procrustes: rmse 0.001043705  max resid 0.008307469 
... Similar to previous best
Run 15 stress 0.1958752 
Run 16 stress 0.1858789 
Run 17 stress 0.1835369 
... Procrustes: rmse 0.001901863  max resid 0.01181043 
Run 18 stress 0.1873897 
Run 19 stress 0.188248 
Run 20 stress 0.1935297 
Run 21 stress 0.1962298 
Run 22 stress 0.1835252 
... Procrustes: rmse 0.000572239  max resid 0.002508577 
... Similar to previous best
Run 23 stress 0.1853343 
Run 24 stress 0.1835247 
... New best solution
... Procrustes: rmse 0.0003363192  max resid 0.002640972 
... Similar to previous best
Run 25 stress 0.1835397 
... Procrustes: rmse 0.001640891  max resid 0.01062063 
Run 26 stress 0.1835251 
... Procrustes: rmse 0.0003441826  max resid 0.002335534 
... Similar to previous best
Run 27 stress 0.1835247 
... Procrustes: rmse 0.0005054195  max resid 0.003808553 
... Similar to previous best
Run 28 stress 0.1851553 
Run 29 stress 0.1853071 
Run 30 stress 0.1881066 
Run 31 stress 0.1885172 
Run 32 stress 0.1835275 
... Procrustes: rmse 0.001084588  max resid 0.008087213 
... Similar to previous best
Run 33 stress 0.1907827 
Run 34 stress 0.1874641 
Run 35 stress 0.1850632 
Run 36 stress 0.1835257 
... Procrustes: rmse 0.0008646151  max resid 0.005770151 
... Similar to previous best
Run 37 stress 0.1835266 
... Procrustes: rmse 0.0008395579  max resid 0.006725346 
... Similar to previous best
Run 38 stress 0.1835264 
... Procrustes: rmse 0.001070495  max resid 0.008509586 
... Similar to previous best
Run 39 stress 0.188025 
Run 40 stress 0.1835247 
... Procrustes: rmse 0.0004901265  max resid 0.003618958 
... Similar to previous best
*** Best solution repeated 8 times

stressplot(df1s.mds3, main="Shepard plot")

df1s.mds3

Call:
metaMDS(comm = df1s.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.1835247 
Stress type 1, weak ties
Best solution was repeated 8 times in 40 tries
The best solution was from try 24 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 2d
df1s.mds2 <- metaMDS(df1s.bc,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.2620287 
Run 1 stress 0.2553686 
... New best solution
... Procrustes: rmse 0.0838107  max resid 0.2759917 
Run 2 stress 0.2598275 
Run 3 stress 0.2655848 
Run 4 stress 0.2633175 
Run 5 stress 0.2697909 
Run 6 stress 0.264019 
Run 7 stress 0.2941411 
Run 8 stress 0.2666523 
Run 9 stress 0.263381 
Run 10 stress 0.2628465 
Run 11 stress 0.2682624 
Run 12 stress 0.2617334 
Run 13 stress 0.2674227 
Run 14 stress 0.2624441 
Run 15 stress 0.2628829 
Run 16 stress 0.2568964 
Run 17 stress 0.266174 
Run 18 stress 0.2700445 
Run 19 stress 0.2595927 
Run 20 stress 0.270985 
Run 21 stress 0.259557 
Run 22 stress 0.2607448 
Run 23 stress 0.2652935 
Run 24 stress 0.2598874 
Run 25 stress 0.271551 
Run 26 stress 0.267075 
Run 27 stress 0.2667411 
Run 28 stress 0.2636827 
Run 29 stress 0.2702476 
Run 30 stress 0.2628844 
Run 31 stress 0.2633379 
Run 32 stress 0.2653752 
Run 33 stress 0.259173 
Run 34 stress 0.2560561 
Run 35 stress 0.258853 
Run 36 stress 0.2587636 
Run 37 stress 0.2872915 
Run 38 stress 0.2617806 
Run 39 stress 0.2559239 
Run 40 stress 0.2700518 
Run 41 stress 0.2667658 
Run 42 stress 0.260049 
Run 43 stress 0.2645332 
Run 44 stress 0.2595304 
Run 45 stress 0.2628461 
Run 46 stress 0.2648597 
Run 47 stress 0.2684644 
Run 48 stress 0.2692849 
Run 49 stress 0.2659958 
Run 50 stress 0.2629114 
Run 51 stress 0.2699549 
Run 52 stress 0.269873 
Run 53 stress 0.255736 
... Procrustes: rmse 0.01997224  max resid 0.1061294 
Run 54 stress 0.2557101 
... Procrustes: rmse 0.01306139  max resid 0.08121592 
Run 55 stress 0.2553434 
... New best solution
... Procrustes: rmse 0.01012858  max resid 0.07959973 
Run 56 stress 0.2622363 
Run 57 stress 0.2735116 
Run 58 stress 0.2690758 
Run 59 stress 0.272355 
Run 60 stress 0.2619185 
Run 61 stress 0.2621461 
Run 62 stress 0.2645114 
Run 63 stress 0.2642568 
Run 64 stress 0.2642039 
Run 65 stress 0.2665846 
Run 66 stress 0.2643829 
Run 67 stress 0.2640766 
Run 68 stress 0.2592044 
Run 69 stress 0.2595033 
Run 70 stress 0.262783 
Run 71 stress 0.2670372 
Run 72 stress 0.2662417 
Run 73 stress 0.2742564 
Run 74 stress 0.2624512 
Run 75 stress 0.2670664 
Run 76 stress 0.262287 
Run 77 stress 0.2616465 
Run 78 stress 0.2620381 
Run 79 stress 0.2687522 
Run 80 stress 0.265521 
*** Best solution was not repeated -- monoMDS stopping criteria:
     9: no. of iterations >= maxit
    71: stress ratio > sratmax

stressplot(df1s.mds2, main="Shepard plot")

df1s.mds2

Call:
metaMDS(comm = df1s.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.2553434 
Stress type 1, weak ties
Best solution was not repeated after 80 tries
The best solution was from try 55 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Stress values suggest 3d plot is a more accurate representation

Need to plot dim 1 vs dim 2 and dim 1 vs dim 3 with tsf groups identified by symbols

Generate graphs

This is a variant of code used in worked example for Ch 16. It does use the graph settings defined in appearance.R, which is called earlier.

a<-as.data.frame(scores(df1s.mds3))
a<-cbind(df[c(1:2)],a)   #Add site names & symbols from original data file
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=tsf, group = tsf) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="TSF",
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = tsf, group = tsf ) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="TSF",
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

How would you fit a linear model to assess fire effects?

Take a look at multivariate dispersion first

df.disp <- betadisper(df1s.bc,df$tsf)
anova(df.disp)

Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value    Pr(>F)    
Groups     2 0.48563 0.242815  32.165 7.573e-11 ***
Residuals 76 0.57374 0.007549                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Suggestion of substantial differences in dispersion. Long TSF sites seem quite similar, those with shorter times are more variable and similar to each other.

Do permanova including tsf

df1s.ado <- adonis2(df1s.bc~tsf,data=df,permutations=999)
print(df1s.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ tsf, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
tsf       2   2.6433 0.11039 4.7151  0.001 ***
Residual 76  21.3028 0.88961                  
Total    78  23.9460 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permanova suggests strong signal from tsf

SIMPER analysis for tsf group differences

df.sim <- simper(df1s[,-(1:2)], df$tsf, permutations=1000)
summary(df.sim)

Contrast: long_short 

     average      sd   ratio     ava     avb cumsum      p    
lgui 0.13247 0.14056 0.94250 0.15918 0.26922  0.172 0.3307    
eul  0.09332 0.12167 0.76700 0.02379 0.16793  0.293 0.0509 .  
dcor 0.08201 0.11054 0.74190 0.06582 0.10672  0.399 0.1129    
amac 0.06242 0.05874 1.06250 0.07712 0.07023  0.480 0.9790    
tnig 0.05948 0.07551 0.78760 0.03804 0.07852  0.557 0.8531    
ldel 0.05168 0.04157 1.24340 0.08687 0.00000  0.624 0.0010 ***
pent 0.05090 0.04718 1.07890 0.10975 0.03392  0.690 0.7093    
htal 0.03859 0.05635 0.68490 0.04599 0.02950  0.740 0.3856    
lwhi 0.03208 0.03938 0.81470 0.04549 0.01261  0.782 0.0150 *  
pspe 0.02827 0.03740 0.75580 0.04834 0.00680  0.819 0.0629 .  
amur 0.02647 0.05251 0.50400 0.04512 0.00000  0.853 0.9111    
rdie 0.02561 0.09803 0.26120 0.00512 0.04275  0.886 0.5025    
ecun 0.02190 0.04456 0.49150 0.03603 0.00000  0.915 0.0430 *  
esax 0.01860 0.04544 0.40940 0.01674 0.01687  0.939 0.1289    
ptex 0.01585 0.04231 0.37450 0.01170 0.01687  0.959 0.2208    
aram 0.01495 0.03676 0.40680 0.01517 0.01261  0.979 0.5724    
vros 0.01279 0.03799 0.33650 0.02216 0.00000  0.995 0.8362    
ppor 0.00363 0.01891 0.19220 0.00620 0.00000  1.000 0.6973    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Contrast: long_medium 

     average      sd   ratio     ava     avb cumsum      p    
lgui 0.10221 0.10333 0.98920 0.15918 0.21705  0.133 0.9710    
amac 0.10139 0.11774 0.86110 0.07712 0.16122  0.264 0.1918    
tnig 0.07029 0.10274 0.68410 0.03804 0.10788  0.356 0.6474    
pent 0.06601 0.06862 0.96190 0.10975 0.08112  0.442 0.0220 *  
amur 0.05952 0.11624 0.51210 0.04512 0.07283  0.519 0.0849 .  
dcor 0.05521 0.09128 0.60480 0.06582 0.03801  0.591 0.8202    
ldel 0.04942 0.03935 1.25580 0.08687 0.00185  0.655 0.0010 ***
eul  0.04752 0.07141 0.66540 0.02379 0.07484  0.716 0.9960    
htal 0.03907 0.06803 0.57430 0.04599 0.02956  0.767 0.3636    
vros 0.03267 0.09018 0.36230 0.02216 0.03599  0.809 0.2148    
lwhi 0.02992 0.03771 0.79360 0.04549 0.00840  0.848 0.0020 ** 
pspe 0.02952 0.03909 0.75530 0.04834 0.00806  0.887 0.0020 ** 
ecun 0.02116 0.04260 0.49670 0.03603 0.00000  0.914 0.0030 ** 
rdie 0.01844 0.06391 0.28850 0.00512 0.02970  0.938 0.8312    
aram 0.01711 0.04409 0.38800 0.01517 0.01739  0.960 0.5035    
esax 0.01138 0.02737 0.41600 0.01674 0.00420  0.975 0.6783    
ptex 0.00983 0.02516 0.39070 0.01170 0.00658  0.988 0.7293    
ppor 0.00929 0.03741 0.24840 0.00620 0.01107  1.000 0.3497    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Contrast: short_medium 

     average      sd   ratio     ava     avb cumsum      p   
lgui 0.15998 0.16074 0.99530 0.26922 0.21705  0.205 0.0100 **
amac 0.10618 0.14346 0.74010 0.07023 0.16122  0.340 0.1489   
eul  0.10520 0.12391 0.84900 0.16793 0.07484  0.475 0.0040 **
tnig 0.08903 0.11469 0.77630 0.07852 0.10788  0.589 0.0709 . 
dcor 0.07378 0.14450 0.51060 0.10672 0.03801  0.683 0.2507   
pent 0.04966 0.08217 0.60440 0.03392 0.08112  0.747 0.8022   
amur 0.04038 0.12310 0.32800 0.00000 0.07283  0.798 0.6943   
rdie 0.03751 0.11223 0.33420 0.04275 0.02970  0.847 0.1548   
htal 0.03122 0.07878 0.39630 0.02950 0.02956  0.886 0.7133   
vros 0.02238 0.09233 0.24240 0.00000 0.03599  0.915 0.6024   
aram 0.01559 0.04775 0.32650 0.01261 0.01739  0.935 0.5914   
ptex 0.01278 0.04462 0.28650 0.01687 0.00658  0.951 0.4216   
esax 0.01161 0.04217 0.27530 0.01687 0.00420  0.966 0.5874   
lwhi 0.01095 0.03778 0.28970 0.01261 0.00840  0.980 0.9930   
pspe 0.00840 0.02396 0.35050 0.00680 0.00806  0.991 0.9990   
ppor 0.00608 0.03495 0.17380 0.00000 0.01107  0.999 0.5564   
ldel 0.00101 0.00583 0.17380 0.00000 0.00185  1.000 1.0000   
ecun 0.00000 0.00000     NaN 0.00000 0.00000  1.000 0.9990   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 1000

Suppose you think of a reptile assemblage as simply the species that are present, regardless of how common they are. What would you conclude about the infuence of time since fire now?

Repeat the preceding analysis after creating a new data file that is presence-absence only.

df1.jac <- vegdist(df1,binary=TRUE,method='jaccard')
#Run 3d first
df1s.mds3 <- metaMDS(df1.jac,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.1404128 
Run 1 stress 0.1478658 
Run 2 stress 0.1416978 
Run 3 stress 0.1405588 
... Procrustes: rmse 0.01500187  max resid 0.08886188 
Run 4 stress 0.1404109 
... New best solution
... Procrustes: rmse 0.002860922  max resid 0.02265594 
Run 5 stress 0.145486 
Run 6 stress 0.1476233 
Run 7 stress 0.1418068 
Run 8 stress 0.1423797 
Run 9 stress 0.1457509 
Run 10 stress 0.1416979 
Run 11 stress 0.1404126 
... Procrustes: rmse 0.002853227  max resid 0.02281661 
Run 12 stress 0.140558 
... Procrustes: rmse 0.01456468  max resid 0.08805689 
Run 13 stress 0.1414712 
Run 14 stress 0.1411415 
Run 15 stress 0.1405617 
... Procrustes: rmse 0.0149327  max resid 0.0899859 
Run 16 stress 0.1415481 
Run 17 stress 0.1404251 
... Procrustes: rmse 0.00228483  max resid 0.01842417 
Run 18 stress 0.1428212 
Run 19 stress 0.1405586 
... Procrustes: rmse 0.01462785  max resid 0.08844371 
Run 20 stress 0.1404136 
... Procrustes: rmse 0.003032611  max resid 0.02427334 
Run 21 stress 0.1405605 
... Procrustes: rmse 0.01485076  max resid 0.08965517 
Run 22 stress 0.1416429 
Run 23 stress 0.1499082 
Run 24 stress 0.1420545 
Run 25 stress 0.1419425 
Run 26 stress 0.1405616 
... Procrustes: rmse 0.01492593  max resid 0.08996847 
Run 27 stress 0.1405633 
... Procrustes: rmse 0.0149849  max resid 0.09014188 
Run 28 stress 0.1404112 
... Procrustes: rmse 0.0001823383  max resid 0.0008661697 
... Similar to previous best
Run 29 stress 0.1418758 
Run 30 stress 0.1485354 
Run 31 stress 0.1419033 
Run 32 stress 0.1458155 
Run 33 stress 0.1404104 
... New best solution
... Procrustes: rmse 0.0002170324  max resid 0.001584222 
... Similar to previous best
Run 34 stress 0.1404102 
... New best solution
... Procrustes: rmse 0.0002560659  max resid 0.001758468 
... Similar to previous best
Run 35 stress 0.1404098 
... New best solution
... Procrustes: rmse 0.001533377  max resid 0.01214581 
Run 36 stress 0.1414897 
Run 37 stress 0.1404105 
... Procrustes: rmse 0.0003876278  max resid 0.003003357 
... Similar to previous best
Run 38 stress 0.1443789 
Run 39 stress 0.1404132 
... Procrustes: rmse 0.002460885  max resid 0.01976426 
Run 40 stress 0.1488112 
*** Best solution repeated 1 times

stressplot(df1s.mds3, main="Shepard plot")

df1s.mds3

Call:
metaMDS(comm = df1.jac, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 3 
Stress:     0.1404098 
Stress type 1, weak ties
Best solution was repeated 1 time in 40 tries
The best solution was from try 35 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 2d
df1s.mds2 <- metaMDS(df1.jac,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.2009509 
Run 1 stress 0.2057315 
Run 2 stress 0.2115578 
Run 3 stress 0.2019149 
Run 4 stress 0.2042377 
Run 5 stress 0.2295765 
Run 6 stress 0.2137388 
Run 7 stress 0.2200574 
Run 8 stress 0.2095372 
Run 9 stress 0.2181162 
Run 10 stress 0.2120083 
Run 11 stress 0.2100096 
Run 12 stress 0.2094365 
Run 13 stress 0.2059035 
Run 14 stress 0.2068707 
Run 15 stress 0.2138507 
Run 16 stress 0.2247804 
Run 17 stress 0.2205491 
Run 18 stress 0.2031129 
Run 19 stress 0.2066478 
Run 20 stress 0.2133806 
Run 21 stress 0.2148553 
Run 22 stress 0.2136219 
Run 23 stress 0.2228445 
Run 24 stress 0.2027397 
Run 25 stress 0.214167 
Run 26 stress 0.2118895 
Run 27 stress 0.2011779 
... Procrustes: rmse 0.005744405  max resid 0.04013641 
Run 28 stress 0.2011764 
... Procrustes: rmse 0.005803269  max resid 0.03983499 
Run 29 stress 0.2031113 
Run 30 stress 0.2172545 
Run 31 stress 0.2109358 
Run 32 stress 0.2081906 
Run 33 stress 0.2095477 
Run 34 stress 0.2156386 
Run 35 stress 0.2184839 
Run 36 stress 0.2011783 
... Procrustes: rmse 0.006120401  max resid 0.03870025 
Run 37 stress 0.2133788 
Run 38 stress 0.2175712 
Run 39 stress 0.2171102 
Run 40 stress 0.2155669 
Run 41 stress 0.209539 
Run 42 stress 0.2014972 
Run 43 stress 0.2135397 
Run 44 stress 0.2148463 
Run 45 stress 0.2145609 
Run 46 stress 0.2042383 
Run 47 stress 0.2122397 
Run 48 stress 0.2221256 
Run 49 stress 0.2075648 
Run 50 stress 0.21128 
Run 51 stress 0.2035834 
Run 52 stress 0.2009593 
... Procrustes: rmse 0.004919144  max resid 0.03284558 
Run 53 stress 0.2262974 
Run 54 stress 0.2096971 
Run 55 stress 0.2233517 
Run 56 stress 0.2040341 
Run 57 stress 0.2081925 
Run 58 stress 0.2031114 
Run 59 stress 0.2030044 
Run 60 stress 0.2035828 
Run 61 stress 0.2045719 
Run 62 stress 0.2135396 
Run 63 stress 0.2209364 
Run 64 stress 0.2156286 
Run 65 stress 0.2015097 
Run 66 stress 0.2168916 
Run 67 stress 0.2036082 
Run 68 stress 0.209978 
Run 69 stress 0.202577 
Run 70 stress 0.2074581 
Run 71 stress 0.2119757 
Run 72 stress 0.2038332 
Run 73 stress 0.2021148 
Run 74 stress 0.2150588 
Run 75 stress 0.21096 
Run 76 stress 0.2260586 
Run 77 stress 0.2030607 
Run 78 stress 0.2131804 
Run 79 stress 0.2133794 
Run 80 stress 0.218075 
*** Best solution was not repeated -- monoMDS stopping criteria:
    13: no. of iterations >= maxit
    67: stress ratio > sratmax

stressplot(df1s.mds2, main="Shepard plot")

df1s.mds2

Call:
metaMDS(comm = df1.jac, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 2 
Stress:     0.2009509 
Stress type 1, weak ties
Best solution was not repeated after 80 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Marginal decision here - stress just over 0.2 for k=2, so might be OK. Stress good for k=3, but 2-d MDS usually better for reporting or showing to audience

a<-as.data.frame(scores(df1s.mds3))
a<-cbind(df[c(1:2)],a)   #Add site names & symbols from original data file
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=tsf, ) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="TSF",
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = tsf, ) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="TSF",
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

Unburned sites separate very clearly on MDS1. Pattern clearer than with abundances included

Examine dispersions

df.disp <- betadisper(df1.jac,df$tsf)
anova(df.disp)

Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value   Pr(>F)   
Groups     2 0.18892 0.094460  6.0078 0.003781 **
Residuals 76 1.19495 0.015723                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Still some difference, though not as pronounced when only pres-abs used

Do permanova including tsf

df1.ado <- adonis2(df1.jac~tsf,data=df,permutations=999)
print(df1.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.jac ~ tsf, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
tsf       2   3.7853 0.20678 9.9057  0.001 ***
Residual 76  14.5210 0.79322                  
Total    78  18.3062 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Clear fire signal

Extension activity

Dixon and her colleagues also recorded a range of habitat variables, which we’ve collected for you in dixonenv.csv

dixonenv <- read_csv("../data/dixonenv.csv")

Rows: 79 Columns: 14── Column specification ────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): site, tsf, veg, aspect
dbl (10): lcwd, shrcov, grcov, litcov, litdep, rocks, elev, warm, cold, twi
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(dixonenv, 10)

They recorded Coarse Woody Debris, which was long-transformed, litter cover (litcov), % cover of groundcover (grcov), and % cover of shrubs (shrcov), and rock cover. In the original paper, Dixon et al. were comfortable that these predictors weren’t correlated.

How is reptile assemblage related to time since fire and these five habitat variables

Check continuous predictors

boxplot(dixonenv$lcwd)

boxplot(dixonenv$litcov)

boxplot(dixonenv$grcov)

boxplot(dixonenv$shrcov)

boxplot(dixonenv$rocks)

#transform rocks to log+1
dixonenv$lrocks <- log10(dixonenv$rocks+1)
# check homog of dispersions on double standardized abundances
df1s.disp <- betadisper(df1s.bc,df$tsf)
anova(df1s.disp)

Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value    Pr(>F)    
Groups     2 0.48563 0.242815  32.165 7.573e-11 ***
Residuals 76 0.57374 0.007549                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

do permanova including tsf plus continuous predictors

df1s.ado <- adonis2(df1s.bc~tsf+lcwd+litcov+grcov+shrcov+lrocks,data=dixonenv,permutations=999)
print(df1s.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ tsf + lcwd + litcov + grcov + shrcov + lrocks, data = dixonenv, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
tsf       2   2.6433 0.11039 4.8306  0.001 ***
lcwd      1   0.4865 0.02032 1.7783  0.051 .  
litcov    1   0.4037 0.01686 1.4757  0.145    
grcov     1   0.3281 0.01370 1.1991  0.297    
shrcov    1   0.3291 0.01374 1.2028  0.305    
lrocks    1   0.3300 0.01378 1.2062  0.266    
Residual 71  19.4253 0.81121                  
Total    78  23.9460 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Run analysis for presence-absence as well
df1.jac.ado <- adonis2(df1.jac~tsf+lcwd+litcov+grcov+shrcov+lrocks,data=dixonenv,permutations=999)
print(df1.jac.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.jac ~ tsf + lcwd + litcov + grcov + shrcov + lrocks, data = dixonenv, permutations = 999)
         Df SumOfSqs      R2       F Pr(>F)    
tsf       2   3.7853 0.20678 10.4082  0.001 ***
lcwd      1   0.3802 0.02077  2.0908  0.044 *  
litcov    1   0.3145 0.01718  1.7296  0.086 .  
grcov     1   0.4097 0.02238  2.2530  0.022 *  
shrcov    1   0.2748 0.01501  1.5115  0.133    
lrocks    1   0.2309 0.01262  1.2700  0.239    
Residual 71  12.9108 0.70527                   
Total    78  18.3062 1.00000                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Main effect is tsf; only some question about CWD when we look at species composition and abundance, while for presence-absence data, there are also effects of CWD and groundcover

mvabund

dfmv <- mvabund(df1)
df.mv <- manyglm(dfmv~dixonenv$tsf+dixonenv$lcwd+dixonenv$litcov+dixonenv$grcov+dixonenv$shrcov+dixonenv$lrocks,family="poisson")
plot(df.mv)

# still hetero vars so use -ve binomial
df.mv1 <- manyglm(dfmv~dixonenv$tsf+dixonenv$lcwd+dixonenv$litcov+dixonenv$grcov+dixonenv$shrcov+dixonenv$lrocks,family="negative_binomial")
plot(df.mv1)

anova(df.mv1)

Time elapsed: 0 hr 1 min 37 sec
Analysis of Deviance Table

Model: dfmv ~ dixonenv$tsf + dixonenv$lcwd + dixonenv$litcov + dixonenv$grcov + dixonenv$shrcov + dixonenv$lrocks

Multivariate test:
                Res.Df Df.diff   Dev Pr(>Dev)    
(Intercept)         78                           
dixonenv$tsf        76       2 331.6    0.001 ***
dixonenv$lcwd       75       1  61.5    0.001 ***
dixonenv$litcov     74       1  25.9    0.287    
dixonenv$grcov      73       1  36.4    0.073 .  
dixonenv$shrcov     72       1  21.5    0.601    
dixonenv$lrocks     71       1  58.0    0.003 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Arguments:
 Test statistics calculated assuming uncorrelated response (for faster computation) 
 P-value calculated using 999 iterations via PIT-trap resampling.

B Simple MDS & Permanova

Let’s return to the Hutto and Barrett (2021) example from the Chapter 15 exercises. There, we used similarity-based analyses to explore the relationship between frog assemblages in ponds and the degree of urbanization surrounding. This seems a question that could just as easily be examined using distances or dissimilarities.

Return to the data and assess whether frog assemblages (using the standardized abundance scale or presence-absence) differ with urbanization.

Did your conclusions differ from those obtained using RDA?

df <- read_csv("../data/huttoamph.csv")

Rows: 50 Columns: 14── Column specification ────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): type
dbl (13): Site, ACRE, BAME, BFOW, GCAR, HCIN, HVER, LCAT, LCLA, LPAL, LSPH, PCRU, PFER
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(df,10)

#ponds 9 and 40 had no frogs. Remove for distance analysis
df <- df %>% 
  filter(Site!= 9 & Site !=40)
df

The explanation of the variables is in the previous chapter’s exercises.

df1 <- df[,-(1:2)]   #Stripped back file without labels, leaving only 
df1.bc <- vegdist(df1,'bray')

#Run 2d first
df1.mds2 <- metaMDS(df1.bc,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.196708 
Run 1 stress 0.1874607 
... New best solution
... Procrustes: rmse 0.06607221  max resid 0.2294443 
Run 2 stress 0.1949763 
Run 3 stress 0.1964022 
Run 4 stress 0.1990471 
Run 5 stress 0.1916105 
Run 6 stress 0.2037482 
Run 7 stress 0.1961229 
Run 8 stress 0.1888276 
Run 9 stress 0.1977393 
Run 10 stress 0.2035879 
Run 11 stress 0.2060001 
Run 12 stress 0.1862237 
... New best solution
... Procrustes: rmse 0.1042587  max resid 0.4865567 
Run 13 stress 0.1888252 
Run 14 stress 0.1920872 
Run 15 stress 0.2075707 
Run 16 stress 0.1854506 
... New best solution
... Procrustes: rmse 0.06781964  max resid 0.2195375 
Run 17 stress 0.1974324 
Run 18 stress 0.1853453 
... New best solution
... Procrustes: rmse 0.006310912  max resid 0.03189332 
Run 19 stress 0.203107 
Run 20 stress 0.1888252 
Run 21 stress 0.1875614 
Run 22 stress 0.1862238 
Run 23 stress 0.1973006 
Run 24 stress 0.1853452 
... New best solution
... Procrustes: rmse 0.0001046814  max resid 0.0004231746 
... Similar to previous best
Run 25 stress 0.2077143 
Run 26 stress 0.1967075 
Run 27 stress 0.1862781 
Run 28 stress 0.1933023 
Run 29 stress 0.1877629 
Run 30 stress 0.1862237 
Run 31 stress 0.2165344 
Run 32 stress 0.1862239 
Run 33 stress 0.1929084 
Run 34 stress 0.196071 
Run 35 stress 0.1967906 
Run 36 stress 0.1990753 
Run 37 stress 0.1862236 
Run 38 stress 0.1888204 
Run 39 stress 0.1877626 
Run 40 stress 0.1920872 
*** Best solution repeated 1 times

stressplot(df1.mds2, main="Shepard plot")

df1.mds2

Call:
metaMDS(comm = df1.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.1853452 
Stress type 1, weak ties
Best solution was repeated 1 time in 40 tries
The best solution was from try 24 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 3d
df.mds3 <- metaMDS(df1.bc,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.129187 
Run 1 stress 0.1384284 
Run 2 stress 0.1384287 
Run 3 stress 0.1289617 
... New best solution
... Procrustes: rmse 0.01116102  max resid 0.05648397 
Run 4 stress 0.128962 
... Procrustes: rmse 0.001080951  max resid 0.006138037 
... Similar to previous best
Run 5 stress 0.1380872 
Run 6 stress 0.1291861 
... Procrustes: rmse 0.01054685  max resid 0.05449119 
Run 7 stress 0.1291862 
... Procrustes: rmse 0.01108717  max resid 0.05660768 
Run 8 stress 0.1289618 
... Procrustes: rmse 0.0009977382  max resid 0.005680432 
... Similar to previous best
Run 9 stress 0.1380886 
Run 10 stress 0.1291601 
... Procrustes: rmse 0.009846495  max resid 0.05162542 
Run 11 stress 0.1289617 
... Procrustes: rmse 0.0001180296  max resid 0.000397725 
... Similar to previous best
Run 12 stress 0.1388068 
Run 13 stress 0.1289614 
... New best solution
... Procrustes: rmse 0.0002094611  max resid 0.001134175 
... Similar to previous best
Run 14 stress 0.1291858 
... Procrustes: rmse 0.01061726  max resid 0.05490047 
Run 15 stress 0.138595 
Run 16 stress 0.1384304 
Run 17 stress 0.1291863 
... Procrustes: rmse 0.01062772  max resid 0.0549278 
Run 18 stress 0.1291602 
... Procrustes: rmse 0.009566248  max resid 0.05100092 
Run 19 stress 0.1291864 
... Procrustes: rmse 0.01111886  max resid 0.05672889 
Run 20 stress 0.1383909 
Run 21 stress 0.1380887 
Run 22 stress 0.1384274 
Run 23 stress 0.1384276 
Run 24 stress 0.128962 
... Procrustes: rmse 0.0002702949  max resid 0.00160678 
... Similar to previous best
Run 25 stress 0.1291856 
... Procrustes: rmse 0.01060586  max resid 0.05479499 
Run 26 stress 0.1289617 
... Procrustes: rmse 0.0002416892  max resid 0.001126536 
... Similar to previous best
Run 27 stress 0.1383939 
Run 28 stress 0.1291857 
... Procrustes: rmse 0.01062561  max resid 0.05487515 
Run 29 stress 0.1380879 
Run 30 stress 0.1291609 
... Procrustes: rmse 0.01000971  max resid 0.05217349 
Run 31 stress 0.1291859 
... Procrustes: rmse 0.01092252  max resid 0.05627383 
Run 32 stress 0.1291858 
... Procrustes: rmse 0.01090207  max resid 0.05607242 
Run 33 stress 0.1291861 
... Procrustes: rmse 0.01051196  max resid 0.05451072 
Run 34 stress 0.1289618 
... Procrustes: rmse 0.0002363175  max resid 0.001405223 
... Similar to previous best
Run 35 stress 0.1291857 
... Procrustes: rmse 0.01088096  max resid 0.05598853 
Run 36 stress 0.1289613 
... New best solution
... Procrustes: rmse 0.0001167837  max resid 0.0004408721 
... Similar to previous best
Run 37 stress 0.1314089 
Run 38 stress 0.1384272 
Run 39 stress 0.1289616 
... Procrustes: rmse 0.0002186302  max resid 0.001234739 
... Similar to previous best
Run 40 stress 0.1291857 
... Procrustes: rmse 0.01067024  max resid 0.0551978 
*** Best solution repeated 2 times

stressplot(df.mds3, main="Shepard plot")

df.mds3

Call:
metaMDS(comm = df1.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.1289613 
Stress type 1, weak ties
Best solution was repeated 2 times in 40 tries
The best solution was from try 36 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

2-d is acceptable

a<-as.data.frame(scores(df1.mds2))
a<-cbind(df[c(1:2)],a)   #Add site names & symbols from original data file
ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=type, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="Urbanization",
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )+
  theme_qk() + scale_color_uchicago()

df1.ado <- adonis2(df1.bc~type,data=df,permutations=999)
print(df1.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.bc ~ type, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)  
type      2    0.912 0.07281 1.7669  0.043 *
Residual 45   11.614 0.92719                
Total    47   12.526 1.00000                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permanova detects a difference; MDS plot suggests low urbanization may be different from the other two, though worth looking at dispersion as well, as spread of this type is less than the other two on plot

Repeat analysis with presence-absence

df1.jac <- vegdist(df1,binary=TRUE,method=‘jaccard’)

#Run 2d first
df1.jac <- vegdist(df1,binary=TRUE,method='jaccard')
df1.mds2 <- metaMDS(df1.jac,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.1551612 
Run 1 stress 0.1837013 
Run 2 stress 0.1696919 
Run 3 stress 0.1880348 
Run 4 stress 0.1635149 
Run 5 stress 0.1630923 
Run 6 stress 0.1585799 
Run 7 stress 0.1536788 
... New best solution
... Procrustes: rmse 0.01900199  max resid 0.06508032 
Run 8 stress 0.1934746 
Run 9 stress 0.1538245 
... Procrustes: rmse 0.01100139  max resid 0.05222997 
Run 10 stress 0.1609337 
Run 11 stress 0.1622269 
Run 12 stress 0.1704331 
Run 13 stress 0.1602261 
Run 14 stress 0.1538244 
... Procrustes: rmse 0.01099654  max resid 0.05221544 
Run 15 stress 0.1696919 
Run 16 stress 0.1621113 
Run 17 stress 0.1630125 
Run 18 stress 0.1585799 
Run 19 stress 0.166129 
Run 20 stress 0.1641746 
Run 21 stress 0.1586729 
Run 22 stress 0.1887895 
Run 23 stress 0.1585798 
Run 24 stress 0.1606139 
Run 25 stress 0.163515 
Run 26 stress 0.1739632 
Run 27 stress 0.1709241 
Run 28 stress 0.1772585 
Run 29 stress 0.1850913 
Run 30 stress 0.1630924 
Run 31 stress 0.1721432 
Run 32 stress 0.1608313 
Run 33 stress 0.1808947 
Run 34 stress 0.1602261 
Run 35 stress 0.1866374 
Run 36 stress 0.1555137 
Run 37 stress 0.1691747 
Run 38 stress 0.1952929 
Run 39 stress 0.1538807 
... Procrustes: rmse 0.008809693  max resid 0.04417994 
Run 40 stress 0.1806084 
Run 41 stress 0.1991962 
Run 42 stress 0.2043037 
Run 43 stress 0.1646708 
Run 44 stress 0.1608262 
Run 45 stress 0.1634767 
Run 46 stress 0.1586728 
Run 47 stress 0.1585799 
Run 48 stress 0.1627099 
Run 49 stress 0.1565718 
Run 50 stress 0.1748722 
Run 51 stress 0.1627089 
Run 52 stress 0.1622264 
Run 53 stress 0.1758648 
Run 54 stress 0.1586729 
Run 55 stress 0.1724591 
Run 56 stress 0.1777698 
Run 57 stress 0.1606135 
Run 58 stress 0.1835377 
Run 59 stress 0.1591044 
Run 60 stress 0.1538245 
... Procrustes: rmse 0.01100665  max resid 0.0522056 
Run 61 stress 0.1591019 
Run 62 stress 0.1732141 
Run 63 stress 0.1641745 
Run 64 stress 0.1704183 
Run 65 stress 0.1535586 
... New best solution
... Procrustes: rmse 0.0107227  max resid 0.05414397 
Run 66 stress 0.1695828 
Run 67 stress 0.1586729 
Run 68 stress 0.1609028 
Run 69 stress 0.1609341 
Run 70 stress 0.1781204 
Run 71 stress 0.1551611 
Run 72 stress 0.1585798 
Run 73 stress 0.1627113 
Run 74 stress 0.1930708 
Run 75 stress 0.1721433 
Run 76 stress 0.1673092 
Run 77 stress 0.1538804 
... Procrustes: rmse 0.01706038  max resid 0.08568356 
Run 78 stress 0.1585799 
Run 79 stress 0.1535586 
... Procrustes: rmse 0.0002906816  max resid 0.001349569 
... Similar to previous best
*** Best solution repeated 1 times

stressplot(df1.mds2, main="Shepard plot")

df1.mds2

Call:
metaMDS(comm = df1.jac, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 2 
Stress:     0.1535586 
Stress type 1, weak ties
Best solution was repeated 1 time in 79 tries
The best solution was from try 65 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 3d
df.mds3 <- metaMDS(df1.jac,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.09466813 
Run 1 stress 0.09466822 
... Procrustes: rmse 0.0005907259  max resid 0.002917344 
... Similar to previous best
Run 2 stress 0.09466795 
... New best solution
... Procrustes: rmse 0.000422125  max resid 0.001876703 
... Similar to previous best
Run 3 stress 0.1097137 
Run 4 stress 0.1090623 
Run 5 stress 0.09466803 
... Procrustes: rmse 7.232359e-05  max resid 0.0003468283 
... Similar to previous best
Run 6 stress 0.1005051 
Run 7 stress 0.09466808 
... Procrustes: rmse 0.0001823469  max resid 0.000917663 
... Similar to previous best
Run 8 stress 0.09693885 
Run 9 stress 0.09466795 
... Procrustes: rmse 0.0001688021  max resid 0.0006306265 
... Similar to previous best
Run 10 stress 0.09466813 
... Procrustes: rmse 0.0002394501  max resid 0.001007719 
... Similar to previous best
Run 11 stress 0.09466792 
... New best solution
... Procrustes: rmse 0.000156739  max resid 0.0005673598 
... Similar to previous best
Run 12 stress 0.100505 
Run 13 stress 0.09466824 
... Procrustes: rmse 0.0001553529  max resid 0.0006098347 
... Similar to previous best
Run 14 stress 0.09466817 
... Procrustes: rmse 0.0001279451  max resid 0.0004343751 
... Similar to previous best
Run 15 stress 0.1005054 
Run 16 stress 0.09466821 
... Procrustes: rmse 0.0001518064  max resid 0.0006174398 
... Similar to previous best
Run 17 stress 0.09466807 
... Procrustes: rmse 0.0001705385  max resid 0.00053697 
... Similar to previous best
Run 18 stress 0.10049 
Run 19 stress 0.1004898 
Run 20 stress 0.09466785 
... New best solution
... Procrustes: rmse 0.0003178359  max resid 0.001519645 
... Similar to previous best
Run 21 stress 0.09693835 
Run 22 stress 0.09466791 
... Procrustes: rmse 0.0003074555  max resid 0.001502487 
... Similar to previous best
Run 23 stress 0.1090588 
Run 24 stress 0.1005046 
Run 25 stress 0.09466806 
... Procrustes: rmse 0.0001456666  max resid 0.0008120412 
... Similar to previous best
Run 26 stress 0.09466794 
... Procrustes: rmse 0.000286181  max resid 0.001024811 
... Similar to previous best
Run 27 stress 0.1005051 
Run 28 stress 0.09466829 
... Procrustes: rmse 0.0004617523  max resid 0.002151037 
... Similar to previous best
Run 29 stress 0.09466836 
... Procrustes: rmse 0.0002740026  max resid 0.001533354 
... Similar to previous best
Run 30 stress 0.09466781 
... New best solution
... Procrustes: rmse 0.0001601026  max resid 0.0007993518 
... Similar to previous best
Run 31 stress 0.09693822 
Run 32 stress 0.09693828 
Run 33 stress 0.09466811 
... Procrustes: rmse 0.0003134413  max resid 0.00135763 
... Similar to previous best
Run 34 stress 0.09466819 
... Procrustes: rmse 0.0003169054  max resid 0.001196281 
... Similar to previous best
Run 35 stress 0.1090587 
Run 36 stress 0.1115488 
Run 37 stress 0.1005051 
Run 38 stress 0.10049 
Run 39 stress 0.09466797 
... Procrustes: rmse 0.000226344  max resid 0.0005714952 
... Similar to previous best
Run 40 stress 0.09466804 
... Procrustes: rmse 0.0002594792  max resid 0.00106972 
... Similar to previous best
*** Best solution repeated 5 times

stressplot(df.mds3, main="Shepard plot")

df.mds3

Call:
metaMDS(comm = df1.jac, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 3 
Stress:     0.09466781 
Stress type 1, weak ties
Best solution was repeated 5 times in 40 tries
The best solution was from try 30 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Again, stress OK with 2D, good with 3D. Look at 2D first

a<-as.data.frame(scores(df1.mds2))
a<-cbind(df[c(1:2)],a)   #Add site names & symbols from original data file
ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=type, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="Urbanization",
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )+
  theme_qk() + scale_color_uchicago()

df1.ado <- adonis2(df1.jac~type,data=df,permutations=999)
print(df1.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.jac ~ type, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)
type      2   0.4525 0.04389 1.0328  0.416
Residual 45   9.8587 0.95611              
Total    47  10.3112 1.00000              

No separation seen on Presence-Absence

C

Griffen et al. (2019) examined the oral microbiome of humans, with a focus on the effects of HIV and AntiRetroviral therapy (ART) on this microbiome. They described the microbiome of 341 patients, who fell into one of two categories: HIV^-, and HIV+ with ART. They also matched their samples as far as possible for sex, along with other conditions that can influence oral microbiomes (Candida infection, current smoking). We’ll use their records for these other categories as well. There were more HIV+ than - subjects, but within these groups, approximately equal numbers of two sexes, two Candida infection status, and two smoking categories.

The bacteriome was recorded separately using 16S RNA sequencing, which overed just over 600 taxa.

The data are available from their paper, as two Excel files in the supplementary information. The metadata gives HIV status, and a raft of demographic information. We’ll just work with HIV, sex, Candida, and current smoking. A second sheet has the bacteriome. The code chunk below reads this file (from a Downloads folder - you’ll need to modify the file location). It also screens the file for any bacterial taxa that were not present in this particular study, which reduces the bacterial taxa to 599.

library(readxl)
#Read in metadata. Lots of information we don't need, so a couple of iterations to get down to a few columns we want - HIV status, candida, current smoking, gender
df <- read_excel("../data/41598_2019_55703_MOESM3_ESM.xlsx", 
     range = "A1:S342", na = "NA")
df<-select(df, HIV, candida, gender, smoking_current)
head(df)

#df1 is the bacterial data
df1 <- read_excel("../data/41598_2019_55703_MOESM2_ESM.xlsx")

New names:

head(df1,10)

df1 <- df1 %>%
  select(where( ~ is.numeric(.x) && sum(.x) != 0))  #Drops any bacterial taxon not recorded in this study (i.e. where it's all zeroes)
df1 <- df1[,-1] #remove col1, which is sample ID

Use a dissimilarity based approach to assess the combined effects of HIV-ART, sex, Candida infection, and current smoking on the bacteriome

Do these factors act independently on the bacteriome?

Which effects are largest (use some graphical methods)?

First steps: Think about standardization and which distance measure you’ll use. Bray-Curtis seems common for these kind of data, and and counts of different taxa vary widely.

df1s <- wisconsin(df1)   #Common to do some standardisation, and counts of different taxa vary widely
df1s.bc <- vegdist(df1s,'bray')   #With standardisation (wisconsin)

Run factorial model using permanova, based on B-C

df1.ado <- adonis2(df1s.bc~HIV*candida*gender*smoking_current,data=df,permutations=999)
print(df1.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ HIV * candida * gender * smoking_current, data = df, permutations = 999)
                                    Df SumOfSqs      R2      F Pr(>F)    
HIV                                  1    1.369 0.01250 4.3776  0.001 ***
candida                              1    1.461 0.01333 4.6694  0.001 ***
gender                               1    0.492 0.00449 1.5737  0.013 *  
smoking_current                      1    1.111 0.01014 3.5513  0.001 ***
HIV:candida                          1    0.318 0.00290 1.0151  0.386    
HIV:gender                           1    0.354 0.00324 1.1332  0.209    
candida:gender                       1    0.257 0.00234 0.8202  0.843    
HIV:smoking_current                  1    0.269 0.00245 0.8584  0.769    
candida:smoking_current              1    0.344 0.00314 1.0989  0.238    
gender:smoking_current               1    0.394 0.00360 1.2602  0.104    
HIV:candida:gender                   1    0.349 0.00319 1.1164  0.234    
HIV:candida:smoking_current          1    0.280 0.00255 0.8943  0.687    
HIV:gender:smoking_current           1    0.265 0.00242 0.8463  0.793    
candida:gender:smoking_current       1    0.254 0.00232 0.8110  0.857    
HIV:candida:gender:smoking_current   1    0.362 0.00330 1.1562  0.190    
Residual                           325  101.664 0.92809                  
Total                              340  109.541 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Model shows main effects, but no interactions important (or at least, detected)

Could you fit a simpler model to the data?

Run simpler model

df2.ado <- adonis2(df1s.bc~HIV+candida+gender+smoking_current,data=df,permutations=999)
print(df2.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ HIV + candida + gender + smoking_current, data = df, permutations = 999)
                 Df SumOfSqs      R2      F Pr(>F)    
HIV               1    1.369 0.01250 4.3775  0.001 ***
candida           1    1.461 0.01333 4.6693  0.001 ***
gender            1    0.492 0.00449 1.5736  0.026 *  
smoking_current   1    1.111 0.01014 3.5512  0.001 ***
Residual        336  105.108 0.95953                  
Total           340  109.541 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

No change; not too surprising, as all four main effects were detected easily.

Now go to data visualization

#Run 3d first
df1s.mds3 <- metaMDS(df1s.bc,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.1573956 
Run 1 stress 0.1630191 
Run 2 stress 0.1576127 
... Procrustes: rmse 0.0118972  max resid 0.1519299 
Run 3 stress 0.1581654 
Run 4 stress 0.1576599 
... Procrustes: rmse 0.009852819  max resid 0.1527284 
Run 5 stress 0.1574955 
... Procrustes: rmse 0.012444  max resid 0.1517661 
Run 6 stress 0.1587469 
Run 7 stress 0.1575933 
... Procrustes: rmse 0.01163693  max resid 0.1521229 
Run 8 stress 0.1579649 
Run 9 stress 0.1578527 
... Procrustes: rmse 0.006323981  max resid 0.09982249 
Run 10 stress 0.1574957 
... Procrustes: rmse 0.01120014  max resid 0.1523462 
Run 11 stress 0.1577498 
... Procrustes: rmse 0.009888467  max resid 0.1527233 
Run 12 stress 0.157409 
... Procrustes: rmse 0.007115967  max resid 0.09726598 
Run 13 stress 0.1573745 
... New best solution
... Procrustes: rmse 0.01073383  max resid 0.1524203 
Run 14 stress 0.1573637 
... New best solution
... Procrustes: rmse 0.009102594  max resid 0.1079765 
Run 15 stress 0.1578146 
... Procrustes: rmse 0.01085179  max resid 0.153061 
Run 16 stress 0.1579197 
Run 17 stress 0.1579066 
Run 18 stress 0.1575399 
... Procrustes: rmse 0.01237825  max resid 0.1521774 
Run 19 stress 0.1575654 
... Procrustes: rmse 0.00849738  max resid 0.1079133 
Run 20 stress 0.1577392 
... Procrustes: rmse 0.01138932  max resid 0.1522989 
Run 21 stress 0.1577813 
... Procrustes: rmse 0.01136102  max resid 0.1528569 
Run 22 stress 0.159205 
Run 23 stress 0.1579391 
Run 24 stress 0.1591617 
Run 25 stress 0.1577727 
... Procrustes: rmse 0.008387684  max resid 0.0967774 
Run 26 stress 0.1579649 
Run 27 stress 0.1582177 
Run 28 stress 0.1576498 
... Procrustes: rmse 0.01137476  max resid 0.1525317 
Run 29 stress 0.1591301 
Run 30 stress 0.157277 
... New best solution
... Procrustes: rmse 0.008559541  max resid 0.1081297 
Run 31 stress 0.1579428 
Run 32 stress 0.1576009 
... Procrustes: rmse 0.01273404  max resid 0.1520734 
Run 33 stress 0.1576118 
... Procrustes: rmse 0.01247834  max resid 0.1523452 
Run 34 stress 0.1582027 
Run 35 stress 0.158822 
Run 36 stress 0.158023 
Run 37 stress 0.1575531 
... Procrustes: rmse 0.007344253  max resid 0.1083239 
Run 38 stress 0.157274 
... New best solution
... Procrustes: rmse 0.0003543994  max resid 0.003824715 
... Similar to previous best
Run 39 stress 0.1578054 
Run 40 stress 0.1578922 
*** Best solution repeated 1 times

stressplot(df1s.mds3, main="Shepard plot")

df1s.mds3

Call:
metaMDS(comm = df1s.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.157274 
Stress type 1, weak ties
Best solution was repeated 1 time in 40 tries
The best solution was from try 38 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 2d
df1s.mds2 <- metaMDS(df1s.bc,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.2249552 
Run 1 stress 0.2251211 
... Procrustes: rmse 0.009373779  max resid 0.1481459 
Run 2 stress 0.2255671 
Run 3 stress 0.2261504 
Run 4 stress 0.2254572 
Run 5 stress 0.2257456 
Run 6 stress 0.2265326 
Run 7 stress 0.2270347 
Run 8 stress 0.2258506 
Run 9 stress 0.2258621 
Run 10 stress 0.225423 
... Procrustes: rmse 0.01541529  max resid 0.2310901 
Run 11 stress 0.2248945 
... New best solution
... Procrustes: rmse 0.00452352  max resid 0.08220424 
Run 12 stress 0.2257927 
Run 13 stress 0.2249545 
... Procrustes: rmse 0.004482882  max resid 0.0822275 
Run 14 stress 0.2253883 
... Procrustes: rmse 0.01627316  max resid 0.2347665 
Run 15 stress 0.224955 
... Procrustes: rmse 0.004520226  max resid 0.0822353 
Run 16 stress 0.2250835 
... Procrustes: rmse 0.008578052  max resid 0.1524037 
Run 17 stress 0.2275692 
Run 18 stress 0.2254386 
Run 19 stress 0.2254065 
Run 20 stress 0.2253846 
... Procrustes: rmse 0.01579237  max resid 0.2345835 
Run 21 stress 0.2248481 
... New best solution
... Procrustes: rmse 0.006292327  max resid 0.08260263 
Run 22 stress 0.2264731 
Run 23 stress 0.2250843 
... Procrustes: rmse 0.01060527  max resid 0.1508133 
Run 24 stress 0.2251181 
... Procrustes: rmse 0.008233582  max resid 0.1483655 
Run 25 stress 0.2255503 
Run 26 stress 0.2261628 
Run 27 stress 0.2251126 
... Procrustes: rmse 0.008253967  max resid 0.1490181 
Run 28 stress 0.2265968 
Run 29 stress 0.2248932 
... Procrustes: rmse 0.006309882  max resid 0.08268836 
Run 30 stress 0.2253832 
Run 31 stress 0.2277665 
Run 32 stress 0.2258463 
Run 33 stress 0.2253349 
... Procrustes: rmse 0.009301725  max resid 0.1493075 
Run 34 stress 0.2256928 
Run 35 stress 0.2254409 
Run 36 stress 0.2248496 
... Procrustes: rmse 0.0002922941  max resid 0.005119062 
... Similar to previous best
Run 37 stress 0.2248936 
... Procrustes: rmse 0.006317961  max resid 0.08246939 
Run 38 stress 0.2260217 
Run 39 stress 0.2273473 
Run 40 stress 0.2255819 
*** Best solution repeated 1 times

stressplot(df1s.mds2, main="Shepard plot")

df1s.mds2

Call:
metaMDS(comm = df1s.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.2248481 
Stress type 1, weak ties
Best solution was repeated 1 time in 40 tries
The best solution was from try 21 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

3D model fits better than 2D. Stress for 2D above 0.2

Start with plots where symbol colour indicates levels of one of the factors

a<-as.data.frame(scores(df1s.mds3))
a<-cbind(df[c(1:4)],a)   #Add site names & symbols from original data file
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=HIV) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = HIV) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="HIV status") & theme_qk() & scale_color_uchicago()

p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=candida) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = candida) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="Candida") & theme_qk() & scale_color_uchicago()

p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=gender) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = gender) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="Sex") & theme_qk() & scale_color_uchicago()

p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=smoking_current) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = smoking_current) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="Current smoking") & theme_qk() & scale_color_uchicago()

Try some visualizations of pairs of factors using MDS1 & MDS2 (while 3D is preferred, stress of 3D is around .15, vs 0.22, so not huge difference) Separate HIV+/-, then use symbol colour to separate second factor.

p1<-ggplot(data=subset(a, candida=="Yes"), aes(x=NMDS1, y=NMDS2, color=HIV) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title="Candida")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=subset(a, candida =="No"), aes(x=NMDS1, y=NMDS2, color = HIV) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title = "No Candida")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )
p_combined <- p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

p_ranges_x <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_x[[1]]$range$range,
  ggplot_build(p_combined[[2]])$layout$panel_scales_x[[1]]$range$range)

p_ranges_y <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_y[[1]]$range$range,
                ggplot_build(p_combined[[2]])$layout$panel_scales_y[[1]]$range$range)
p_combined & 
  xlim(min(p_ranges_x), max(p_ranges_x)) & 
  ylim(min(p_ranges_y), max(p_ranges_y))

p1<-ggplot(data=subset(a, HIV=="Yes"), aes(x=NMDS1, y=NMDS2, color=candida) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title="HIV+")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=subset(a, HIV =="No"), aes(x=NMDS1, y=NMDS2, color = candida) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title = "HIV-")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )
p_combined <- p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

p_ranges_x <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_x[[1]]$range$range,
  ggplot_build(p_combined[[2]])$layout$panel_scales_x[[1]]$range$range)

p_ranges_y <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_y[[1]]$range$range,
                ggplot_build(p_combined[[2]])$layout$panel_scales_y[[1]]$range$range)
p_combined & 
  xlim(min(p_ranges_x), max(p_ranges_x)) & 
  ylim(min(p_ranges_y), max(p_ranges_y))

p1<-ggplot(data=subset(a, HIV=="Yes"), aes(x=NMDS1, y=NMDS2, color = smoking_current) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title="HIV+")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=subset(a, HIV =="No"), aes(x=NMDS1, y=NMDS2, color = smoking_current) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title = "HIV-")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )
p_combined <- p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

p_ranges_x <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_x[[1]]$range$range,
  ggplot_build(p_combined[[2]])$layout$panel_scales_x[[1]]$range$range)

p_ranges_y <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_y[[1]]$range$range,
                ggplot_build(p_combined[[2]])$layout$panel_scales_y[[1]]$range$range)
p_combined & 
  xlim(min(p_ranges_x), max(p_ranges_x)) & 
  ylim(min(p_ranges_y), max(p_ranges_y))

References

Dixon, Kelly M., Geoffrey J. Cary, Graeme L. Worboys, and Philip Gibbons. 2018. “The Disproportionate Importance of Long‐unburned Forests and Woodlands for Reptiles.” Ecology and Evolution 8 (22): 10952–63. https://doi.org/gfs587.

Griffen, Ann L., Zachary A. Thompson, Clifford J. Beall, Elizabeth A. Lilly, Carolina Granada, Kelly D. Treas, Kenneth R. DuBois, et al. 2019. “Significant Effect of HIV/HAART on Oral Microbiota Using Multivariate Analysis.” Scientific Reports 9 (1): 19946. https://doi.org/gsfpz3.

Hutto, David, and Kyle Barrett. 2021. “Do Urban Open Spaces Provide Refugia for Frogs in Urban Environments?” Plos One 16 (1). https://doi.org/gr2r9n.